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Abstract 

Zitterbewegung, a force-free trembling motion first predicted for relativistic fermions like elec- 
trons [I], was an unexpected consequence of the Dirac equation's unification of quantum mechanics 
and special relativity. Though the oscillatory motion's large frequency and small amplitude have 
precluded its measurement with electrons, zitterbewegung is observable via quantum simulation. 
We engineered an environment for ^'^Rb Bose-Einstein condensates where the constituent atoms 
behaved like relativistic particles subject to the one-dimensional Dirac equation. With direct imag- 
ing, we observed the sub-micrometer trembling motion of these clouds, demonstrating the utility 
of neutral ultracold quantum gases for simulating Dirac particles. 
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Among the great discoveries of the Enhghtenment was the reahzation that physical laws 
are equivalent in all places, at all times, and for all scales; this remains a central tenet in 
contemporary science. Quantum simulation exploits this universality to study the behavior 
of systems that are difficult to access or impossible to manipulate, by performing direct 
measurements on analogue systems composed of well-characterized and highly manipulable 
quantum building blocks. In this work, we used neutral rubidium atoms to simulate zitterbe- 
wegung, a trembling motion usually associated with relativistic electrons, and we illuminate 
its microscopic origins by drawing an analogy to the well-understood atomic physics of Rabi 
oscillations. The Dirac equation - describing the motion of free fermions - is an essential 
part of our current description of nature; by engineering new Dirac particles in novel set- 
tings, we expose the equation's properties by direct measurement. Simulations of the Dirac 
equation have been proposed for superconductors [2] , semiconductors [3l E] , graphene [5] , 
cold atoms [U I614T3] , and photonic systems |13] ; and have been realized with cold atoms [I5] , 
trapped ions (TB] , and photons [T7] . The ion and photon experiments demonstrated zitterbe- 
wegung for quantities analogous to position or time in the Dirac equation. Here, we directly 
observed a neutral-atom BEC undergoing zitterbewegung in space and time. 

The precise control and direct measurement techniques available in systems of ultra- 
cold atoms, coupled with their accessible length and energy scales, make these systems 
ideal for quantum simulation. In this experiment, our quantum building blocks were Bose- 
condensed ^^Rb atoms. Using two counter-propagating Raman lasers (Fig. lA) with wave- 
length A = 790.1 nm, we coupled the atoms' \f = l,mp = =1=1) = IfJ,) atomic hyperfine 
states (comprising our effective two level system) to their external motion [THl [IS] with a 
four-photon Raman transition (Fig. IB). In this environment, each atom's behavior was gov- 
erned by the one-dimensional Dirac Hamiltonian, making its motion analogous to that of a 
relativistic electron. The system's characteristic momentum hk^ = 2'n-h/ X - that of a single 
photon - specifies the recoil energy En = h^k'^/2m = h x 3.68 kHz, where m is the atomic 
mass. These recoil units set the scale for all physical quantities in our analogue system. The 
Raman lasers drove the four-photon |'|", hk^ = Px + 2/iA;p{) -H- \l, hk^ = Px — 2/iA;r) transition 
(resonant when px = 0), where hkx is the atomic momentum along e^; and px will play the 
role of momentum in the Dirac equation. The simulated speed of light c* = 2v-£i = 11.6 mm/s 
was twice the atoms' recoil velocity, a factor of ~ 10^° less than the true speed of light. The 
artificial rest energy m*c*^ = hfl/2 < 1 was a factor of ~ 10^'' less than the electron's 
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rest energy {HQ is the four-photon laser couphng strength). The effective Compton wave- 
length Aq = h/m*c* ~ 1 |J.m, the approximate amplitude of zitterbewegung, exceeded that 
of an electron by a factor of ~ 10^. These new scales enabled our direct measurement of 
zitterbewegung. 

The dynamics of our ultracold ^^Rb atoms were described by the one-dimensional Dirac 
equation [12] 

H-D = {c*p^az + m*c*^(T^) , (1) 

where is the momentum operator; (jx,y,z are the Pauli spin operators; and \ip) is repre- 
sented as a two-component spinor, whose components are defined by |ti)Px), the m* = 
eigenstates of H-d- For the massless, m* = case, this equation simply describes particles 
(positive energy) or anti-particles (negative energy) traveling with velocity ±c*, respectively, 
as depicted by the dashed lines in Fig. IC. The mass term couples together these m* = 
states, producing an avoided crossing (solid curves in Fig. IC) with energy given by the 
familiar relativistic dispersion E(px) = ±{plc*'^ -|- m*^c*^)^/^, gapped ai px = by twice the 
rest energy. In our atomic analogue, the two massless states coupled by the effective rest 
energy physically corresponded to the atomic states moving with velocity ±c*. 

Zitterbewegung arises because the Pauli matrices associated with the two terms in the 
Dirac equation do not commute. In the Heisenberg representation of quantum mechanics 
the operators, not the wavefunctions, depend on time: for example Vx = dx/dt = [x, Hb]/Ht- 
In this formalism, the velocity operator obeys the differential equation 

d /~\9 2^^C /e-\\ 

+ il Vx= ^ Px^x- (2) 

For an initial state |t, = 0), which gives initial conditions (-O^^) = c* and {dvx/dt) = 0, the 
expectation values of the position and velocity observables oscillate with the zitterbewegung 
frequency Q according to 

{x{t)) =x{0) + ^sm{nt) (3) 
An 

{vx{t)) = c*cosm. (4) 



Initial states with {px) 7^ 0, or localized wave packets, follow more complex trajectories [20] . 

Our experiment began with Bose-condensed atoms moving at c*, in the It? Pa; = 0) state, 
physically the atomic Imp = —^,kx = 2/cr) state. The atoms were confined in a harmonic 
trap [{fx,fy,fz) = (38,38,130) Hz] with characteristic timescales greatly exceeding those of 



the zitterbewegung. We suddenly effected the Dirac Hamiltonian by turning on the Raman 
lasers for a variable evolution time t. We measured the system either by imaging the atoms 
immediately following this evolution (to determine the atoms' position) or by releasing the 
atoms from their trap and simultaneously turning off the Raman lasers, allowing for a short 
time-of-flight (TOF, with duration tTor) before imaging (to determine the atoms' velocity). 
Figure 2 A shows the evolution of the signal for several times of flight, and Figs. 2B-D show 
in situ and TOF (ixoF = 550 \is) measurements at several coupling strengths; the velocity- 
dominated TOF images clearly show the expected cosinusoidal behavior. 

This simple description of zitterbewegung assumes that the range 6px of occupied mo- 
mentum states is small compared to 171*0*, and only those states near the avoided-crossing 
structure are populated. To maintain a sufficiently narrow 6px, the spatial size of the sys- 
tem must be at least 6x ^ A^, which, as observed in Ref. [21], is larger than the A^/47r 
amplitude of the zitterbewegung itself. We satisfied this requirement in our experiment 
by using clouds whose Thomas-Fermi radii = 12(2) |am greatly exceeded the measured 
sub-micron zitterbewegung oscillations, and overcame the fundamental measurement chal- 
lenge with good statistics. Just before initializing zitterbewegung, we measured the initial 
position of the BEG by out-coupling and imaging ^ 5 x 10^ atoms |I9]. In principle, 
this allowed us to measure the center of the distribution with an uncertainty estimated by 
Rx/V5 X 103 ^ 0.17 (xm. Our actual measurements, which include technical noise and are 
averages of four independent images, have a typical 0.3 [im rms uncertainty [much less than 
both the distribution's 12(2) (xm width and our 1.75 [im imaging resolution]. 

From fits to data as in Fig. 2B-D - with parameters joint between each in situ and TOF 
pair - we extracted the frequency Q, amplitude A^/47r, and velocity c* of the observed 
zitterbewegung (shown in Fig. 3, see also [IS])- The observed values are attenuated by 
approximately 2.5 from those predicted by Eqs. ^ and 

For our finite-temperature system, higher momentum states are thermally occupied in the 
initial equilibrium system (schematically illustrated in Fig. IC and observed in Fig. 4). The 
zitterbewegung frequency for these states is increased, and the oscillations correspondingly 
dephase, decreasing the observed amplitude. These finite temperature effects give rise to a 
non-participating fraction fk of the atomic population drifting at c*. Indeed, Fig. 4 shows 
that the majority of the "thermal" population surrounding the initial BEG is unaffected by 
the coupling. Additionally, owing to imperfect preparation of the initial \t,Px = 0) state, 
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a fraction /o = 0.15(10) remained at rest [19]. Fluctuations in the background magnetic 
field also contribute to both /o and fk- The rest fraction /q was determined from long TOF 
images such as those in Fig. 4, and the drifting fraction was found by fitting a model 
which included an overall slope (removed in the data of Fig. 2) that specifies fk- The average 
of the extracted fk (shown in Fig. 3B) was used in the theory curve in Fig. 3C to show the 
expected amphtude of in situ oscillations. 

By engineering a two-level quantum system from initial states with opposite velocity, 
we reinterpret the "curious" physics of zitterbewegung in analogy to the Rabi oscillations 
ubiquitous in atomic physics. In this language, the particles trembled because the initial 
state was not an eigenstate of the coupled system; once subject to the Dirac Hamiltonian, 
the system Rabi-oscillated between bare states of equal and opposite velocity. As the atoms' 
coupling was provided by resonant laser light instead of the electrons' rest energy, it is natural 
to think of a Rabi oscillation picture where the mass (coupling) is suddenly turned on and 
off. (Somewhat amusingly, the mechanism by which our laser field (a coherent state of light) 
generates mass is analogous to the Higgs mechanism where a Higgs condensate (a coherent 
matter wave) generates mass in the standard model [22] •) The zitterbewegung of electrons 
arises because two states - particle and antiparticle states - are coupled, and the resulting 
eigenstates are superpositions of the two. Projections of bare electron states onto this basis 
result, as in the case of the atoms, in oscillations between states of opposite velocity. This 
straightforward analogy compels us to accept that the rest energy acts exactly as a coupling 
field and mixes the particle and antiparticle states into eigenstates that are superpositions 
of the two. 

Zitterbewegung, as observed here, is an example of a broader class of phenomena where a 
group of states with differing velocities are quantum mechanically coupled together [23ti25] . 
As with the present case, eigenstates of the full Hamiltonian are static, but superpositions 
can tremble. Neutrino oscillations [26] are an example of this generalization: neutrinos 
are produced by the weak nuclear force in superpositions of the propagating (i.e., mass) 
eigenstates, each with a different mass and, therefore, velocity. While the Dirac equation 
generally applies only to fermionic systems in nature, quantum simulations such as ours 
directly realize Dirac-boson systems in the laboratory [HI EZ], permitting access to new 
classes of experimental systems. While BECs near these Dirac points are short-lived [2511281 - 
[30] . strong interactions, as are present near the superfiuid-Mott transition in an optical 
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lattice, can stably populate these states [SU |32] , for example leading to bosonic composite- 
fermion states [5^1 
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Fig. 1. A Schematic laser geometry. Two counter propagating laser beams (red and blue) 
coupled the Zeeman levels of the ^''Rb EEC's / = 1 ground state. B Coupled states of ^^Rb. 
The Iti) = \Tnp = =pl) states were laser coupled with a four-photon Raman transition, 
connecting states differing in velocity by Av^. C Top panel: Dirac dispersion relationship 
for massless particles (dashed lines) for massive particles such as electrons and positrons 
(solid curves). With suitable values of m* and c*, this same dispersion relationship and its 
underlying Dirac Hamiltonian equally describe relativistic electrons and our atomic system. 
In the vicinity of the depicted avoided crossing, atoms in |'|") move with velocities near 2fR, 
and those in \D have velocities near —2v^. Bottom panel: Typical momentum distribution 
of the BEC (narrow peak) and thermal cloud (broad) in our system. The vertical axis is 
truncated to show detail - the central peak reaches a value of 18 on this scale. 
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A Signal versus fjop. ^'^jc = 3.0 kHz B in-situ and TOF, a/In = 2.4 kHz 
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Fig. 2 Direct detection of zitterbewegung in position and velocity. A We measured the time- 
evolving position of the cloud subject to the Dirac Hamiltonian for various times of flight. 
We plot in situ (tTOF = 0) and TOF data, showing how the data transform from reflecting 
the atoms' position to indicating their velocity, with f2 = 27r x 3.0 kHz. B-D We repeated 
the in situ (red) and tTOF = 550 [is (black) measurements for various coupling strengths. 
Shown are (B) Q = 2n x 2.4 kHz, (C) = 27r x 3.1 kHz, and (D) Q = 27i x 4.3 kHz. 
Statistical uncertainties are shown with typical error bars in each set of data. 
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A measured frequency B measured fraction, fi^ C measured amplitude with 
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Fig. 3. Observed zitterbewegung parameters A Coupling frequency f2/27r measured for dif- 
ferent values of the Raman laser intensity. As a four-photon transition, this scales quadrat- 
ically with the laser intensity (grey curve). Statistical uncertainties are smaller than the 
symbol size, and the observed scatter likely results from fluctuations in the bias magnetic 
field [T^. B Fraction of atoms fk not participating in zitterbewegung, and their average 
(line). Uncertainties are statistical, and derived from the fit that determines the fraction. 
C Observed in situ amplitude of zitterbewegung (symbols), plotted along with our model, 
including the effects of non-participating atoms with fractions fk and /o = 15(10)%. Uncer- 
tainties in measured values are statistical. 
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Fig. 4. Non-participating atoms. These 30 ms TOF images, which reflect the momentum 
distribution of the atoms, were recorded at three evolution times: A Qt ^ 0, B 7r/4 
and C r2t ~ 7r/2. In C, the transfer of the atoms from thermal cloud was significant only 
near Px = 0, leaving the majority of the thermal atoms behind, thus contributing to the 
non-participating fraction fk- These images show much larger BECs than were used for 
the in situ and short TOF measurements of Figs. 2 and 3, and have a correspondingly 
larger condensate fraction. The atoms in the mp — state result from imperfect initial 
preparation, and contribute to the fraction /q. 
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SUPPLEMENTARY MATERIAL 
System preparation 

Our experiments began with N 5 x 10^ atom optically-trapped ^^Rb BECs [fc — 
0.75(10) condensate fraction] in the |/ = l,mF — —1) ground state, subject to a uniform 
Bq — 2.1 mT bias magnetic field. We transferred these atoms (at rest) to |/ = l,mF — 0) 
using an adiabatic rapid passage technique; a fixed frequency 15.0 MHz radiofrequency 
magnetic field coupled the different uif states together as the bias magnetic field was 
swept through resonance. Using a pair of Raman beams counterpropagating along e^; 
with wavelength A = 790.1 nm and frequency difference Sou — Qf^^bBo + 4£'r + e (where 
e — h X 32 kHz is the quadratic Zeeman shift), a 30 |as 7r-pulse transferred approximately 
85% of the atoms from \mF — 0,kx — 0) to \mF — —i,kx — 2/cr) (moving with velocity 
V = 2hk-^/m = c*). Before the trap appreciably altered their velocity (200 [is), we changed 
the Raman laser's frequency difference to 6u = gF[iBBo, bringing \mF = —1, kx = 2/cr) and 
\mF = +^,kx — —2k-p) into four-photon resonance. We then suddenly introduced a four- 
photon Raman coupling between these states (Fig. SI), and allowed the system to evolve 
under this new Hamiltonian for an evolution time t. 

Imaging 

Just before transferring the EEC into |mp = —l,k,j. = 2k^), two 6.8 GHz microwave 
pulses spaced in time by 50 ms each out-coupled ~ 10% of the atoms to the / = 2 hyperfine 
manifold. These atoms were separately imaged (without repumping on the f — 1 to f — 2 
transition) leaving the atoms in / = 1 undisturbed. These f — 2 atoms served two purposes: 
1) by setting the microwave frequency 2 kHz above (first pulse) and 2 kHz below (second 
pulse) resonance, we tracked shifts in the bias field that would change our four-photon 
Raman resonance condition. Upon analyzing the data, we rejected points where the atom 
number difference between these two images was greater than two standard deviations from 
being equal; 2) we determined the BEC's position immediately before each zitterbewegung 
experiment began, allowing us to cancel shot-to-shot variations in the trap position. The 
beginning of the three transfer pulses - two microwave outcoupling pulses, and the final 
four-photon Raman pulse - were each separated in time by 50 ms. As three periods of a 
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60 Hz cycle, this separation was chosen to reduce magnetic field background fiuctuations at 
the power line frequency, and to facilitate rethermalization between pulses. 

These experiments implemented either in situ or time-of-flight (TOF) absorption imaging. 
For in situ measurements, the Raman and trapping beams remained on during the 40 |J.s 
absorption imaging pulses. For timc-of-fiight, these were removed during TOF during which 
time the atoms flew ballistically for ixoF and were subsequently absorption-imaged. We used 
high intensity imaging, with intensity / ~ 3/sat (where Igat is the saturation intensity), that 
reduced the effective optical depth (SI) and gave better signal-to-noise in the determination 
of the clouds' positions. 



Atomic Dirac Hamiltonian 

Atomic dynamics along e^; of the three uif states comprising the 5Si/2, f — i electronic 
ground state manifold subject to our 2-photon Raman process (Fig. SlA) are described by 
the 3-level Hamiltonian 

" I ^ J — m — ^^'^ ~ (^z,z)ne. (5) 

Eji = h^k'^/2m and hk-^ = 27r^/A are the single photon recoil energy and momentum 
associated with the Raman laser's wavelength A; m is the atomic mass, VL2 is the 2-photon 
Raman coupling strength; e is the quadratic Zeeman shift that energetically displaces the 
mp — ^ state; a^^z are the generalized Pauli operators for a spin-1 system; and I3 is the 
3x3 identity. The eigenenergies of this system are plotted in Fig. SIB. We concentrate on 
the avoided crossing k^ — ^ between the states that adiabatically connect to Im^ = — 1) 
and \mF — +1). By adiabatically eliminating the lowest-energy eigenstate, we obtain the 
effective 2x2 Hamiltonian, 



H 



2x2 



h\kl + ^kl) 
2^ + 



1 H O-x + -TT^^z, (6) 

m 2 



which includes a global rotation of the system Ux ^ c^y, CFy y d'z-i ^z ^ ^x- 

For kx/kK < 1, 

the effective coupling is Q = /j.Q|/2(4£'i — he) and E4 — ElQ,/Q,2- Ignoring the uniform 
energy offset, we identify the parameters from the Dirac Hamiltonian (Eqn. 1): the effective 
c* = 2hk^/m is twice the atomic recoil velocity, the rest energy m*c*^ — Kl/2 is the couphng 
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strength, and the Compton wavelength = h/m*c* = Snhk-n/mQ sets the scale for the 
zitterbewegung's amplitude. 




momentum {ki} 



FIG. 1. Detailed experimental setup A Schematic of laser set up. Two Raman laser beams 
(orange and yellow) counterpropagating along ex couple the three Zeeman levels of ^^Rb's / = 1 
electronic ground state. This coupling connects states differing in momentum by ±2^^^; here, an 
atom initially in Imp = 0,kx) couples to \mp = ^l,kx ± = |t)i)- B Full, 3-level dispersion 
relationship described by Eqn. SI with e = O-Er and ft = 3E^. Red indicates mp = —1, black 
mp = 0, and blue mp = +1. Dashed curves show the uncoupled "bare" states while solid curves 
show the dressed states, with color indicating the dominant contribution to the superposition. We 
concentrate on the high-energy avoided crossing, which is shown in greater detail in the lower left. 
This region of the dispersion relationship simulates the ID Dirac Hamiltonian. 



Analysis 

Our data consisted of density distributions n{x, y), integrated along by the absorption 
imaging process. For each set of experimental parameters, one to five images were recorded 
(more for the smaller- signal small-TOF data). To determine the center-of-mass position 
of each image, the region for which the density was greater than ^ 5% of the peak was 
determined, and the region outside of this was masked (our results are insensitive to the 
exact value of the threshold). The center-of-mass of the density distribution was deter- 
mined by numerically integrating the distribution to obtain x = xn{x , y) / n{x , y) and 
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y = ^ ?/n(x, y)/ ^ n(x, y), where the sums ran over all pixellated positions x,y inside the 
measurement region of each digital image. We made these measurements at various times 
of flight tTOF, and used the data to determine expectation values for both the position and 
velocity of the trapped atoms. 

In our experiments, the less-than-perfect transfer of atoms between internal states con- 
tributes to the reduction of the zitterbewegung amplitude. We use the parameter /o to 
represent the fraction of atoms remaining in the initial state Imp = 0,kx) following trans- 
fer to \mF = —^,kx + 2k-£i). These atoms remain at rest and decrease the overall ampli- 
tude of the signal in position and velocity. The parameter represents the atoms not 
participating in the four-photon process that transfers atoms from Imp = —1, + 2/;;^) to 
\mp = 1, kx — 2k-£i). These atoms moved with speed c* and contribute a background slope 
to the zitterbewegung signal. To account for both types of non-participating atoms, we fit 
the time series to the function 



x{t,tTOF) = (1 - /o) 



sin(fit + 0o) , , , A ^ 
^ \- tTOF cos{ilt + (po) 



+ fkC*t + Xo 



(7) 



where xq is the initial offset position. From independent (long-TOF) measurements, we 
determined /o to be approximately 15%, and we fix this value in the fits. Using two sets of 
data, one in situ and one with txcF = 550 |xs, we performed joint fits for each same laser 
intensity (four-photon coupling strength). In the initial analysis, we fix c* = 2hk^/m and fit 
the data to extract the parameters fk, ^, 0o and Xq. Next, we found the average of fk as a 
function of Q to use in the model. Finally, we remove the background slope due to the fkC*t 
term from the same five sets of data using a simple linear fit. We refit the remaining signal 
to a "fully participating" model [Eq. ([T]) with fo = fk = 0] fixing n (from the original fit) to 
extract the effective speed of light parameter c*. We found the zitterbewegung amplitude 
(of the participating atoms) c*/Q. Uncertainties in these parameters were statistical, and 
were derived by the fitting program from the statistical uncertainties in the in situ and TOF 
data. The model used to predict the amplitude is given by (1 — /o)(l — fk)2hk^/mQ, and 
the uncertainty is dominated by our systematic uncertainty in fo, which is due to magnetic 
field variations. 

The background slope due to the fkC*t term was subtracted from the data presented in 
Fig. 2. The curves are calculated from the values from the original fit using Eq. [7] without 
the fkC*t term. 

15 



Supplementeiry references 

SI. G. Reinaudi, T. Lahaye, Z. Wang, D. Guery-Odelin, Opt. Lett. 32, 3143 (2007). 



16 



